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Satellite observations have shown a positive correlation between cloud amount and 
aerosol optical thickness (AOT) that can be explained by the humidification of aerosols 
near clouds, and/or by cloud contamination by sub-pixel size clouds and the cloud 
adjacency effect. The last effect may substantially increase reflected radiation in cloud-free 
columns, leading to overestimates in the retrieved AOT. For clear-sky areas near boundary 
layer clouds the main contribution to the enhancement of clear sky reflectance at shorter 
wavelengths comes from the radiation scattered into clear areas by clouds and then 
scattered to the sensor by air molecules. Because of the wavelength dependence of air 
molecule scattering, this process leads to a larger reflectance increase at shorter 
wavelengths, and can be corrected using a simple two-layer model [18]. However, 
correcting only for molecular scattering skews spectral properties of the retrieved AOT. 
Kassianov and Ovtchinnikov [9] proposed a technique that uses spectral reflectance ratios 
to retrieve AOT in the vicinity of clouds; they assumed that the cloud adjacency effect 
influences the spectral ratio between reflectances at two wavelengths less than it 
influences the reflectances themselves. This paper combines the two approaches: It 
assumes that the 3D correction for the shortest wavelength is known with some 
uncertainties, and then it estimates the 3D correction for longer wavelengths using a 
modified ratio method. The new approach is tested with 3D radiances simulated for 26 
cumulus fields from Large-Eddy Simulations, supplemented with 40 aerosol profiles. The 
results showed that (i) for a variety of cumulus cloud scenes and aerosol profiles over 
ocean the 3D correction due to cloud adjacency effect can be extended from shorter to 
longer wavelengths and (ii) the 3D corrections for longer wavelengths are not very 
sensitive to unbiased random uncertainties in the 3D corrections at shorter wavelengths. 

Published by Elsevier Ltd. 


1. Introduction 

There is a positive correlation between cloud amount 
and aerosol optical thickness (AOT) (e.g., [16,8,31,19]). 
This correlation can be explained by the humidification 
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of aerosols in the moist cloud environment [23,1 ] and/or 
by a transition between aerosol and clouds where the 
distinction between cloudy and cloud-free air becomes 
problematic [13]. The correlation is also partly due to 
remote sensing artifacts such as cloud contamination 
(e.g., [31,14,27]). 

For low clouds, both sub-pixel cloud contamination 
(e.g., [24,31]) and the cloud adjacency effect (e.g., [12,28,26]) 
substantially increase reflected radiation, thus leading to 
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significant overestimates of the retrieved AOT. In this 
paper we focus on the cloud adjacency effect assuming 
cloud-free columns. 

In the vicinity of boundary layer clouds over ocean, the 
main contribution to the enhancement of clear sky reflec- 
tances at shorter wavelengths comes from the radiation 
scattered into clear areas by clouds and then scattered to 
the sensor by air molecules [29]. This leads to a larger 
reflectance increase for shorter wavelengths, or to a “bluing” 
of aerosols near clouds. Marshak et al. [18] proposed a simple 
two-layer model to correct for the contribution from 
enhanced molecular scattering that results from the pre- 
sence of nearby clouds. This model was recently applied to a 
full Moderate Resolution Imaging Spectroradiometer 
(MODIS) granule [30]. However, correcting only for molecu- 
lar scattering skews spectral properties of the retrieved AOT 
and overcorrects the Angstrom exponent. 

Following a different approach, [9] proposed and 
further developed [10,11] a technique that uses spectral 
reflectance ratios to retrieve AOT in the vicinity of clouds. 
Their main assumption is that the cloud adjacency effect 
impacts the spectral ratio between reflectances at two 
wavelengths less than it impacts the reflectances them- 
selves. Thus, the ratio of 3D radiances at two wavelengths 
X\ and A 2 is approximately equal to the ratio of their ID 
counterparts, i.e. 

^3d(M )/ ^30(^2) ~ KidUi )/ ^10(^2) = p(h > h) ( 1 ) 

where R 3D (A) and R 1D (A) are respectively 3D and ID 
radiances at wavelength A. With three wavelengths (470, 
660 and 870 nm), the ratio method uses two independent 
ratios from which the two parameters of the AOT spectral 
dependence described by a power-law are retrieved. 

Here we combine the two approaches. We assume that 
the 3D correction for the shortest wavelength (466 nm in 
this study) is known (with some level of uncertainty). 
Based on it, we estimate the 3D correction for longer 
wavelengths (855 nm in this study) using a modified ratio 
method. We test our new approach of extending the 
correction from shorter to longer wavelengths with 3D 
(and ID) radiances calculated by the Spherical Harmonic 
Discrete Ordinate Method (SHDOM) for atmospheric radia- 
tive transfer [7] for 26 cumulus fields from the UCLA 
Large-Eddy Simulations (LES) model [22], supplemented 
with 40 aerosol profiles based on GEOS-5 [21] global 
reanalysis. Sections 2 and 3 describe our approach and 
dataset, respectively, Section 4 presents the results, and 
Section 5 offers a summary and discussion. 

2. Approach 

Let the observed (normalized) radiance R 3 D (A) at wave- 
length A be represented as a sum of two components: 
Rid(A) and 4(A), i.e. 

R3d( 2) = ^io(2)+4(A). (2a) 

here Ri D (A) is a ID counterpart of R 3D (A) that accounts only 
for a clear-sky column and ignores any cloud-related 
adjacency effects. The difference between R 3D (A) and 
Rid(A) is 4(A) and will be called a spectral 3D correction. 
The question we address here is how the 3D correction for 


longer wavelengths can be estimated from the shorter 
wavelength correction. For example, if 4(A!=466nm) is 
known, how 4(A 2 =855 nm) can be assessed. The impor- 
tance of using shorter wavelength correction for estimat- 
ing its longer wavelength counterparts follows from our 
ability to account for the cloud-Rayleigh radiation interac- 
tion [18] that dominates for shorter wavelengths over dark 
surfaces with an aerosol layer predominantly below the 
cloud tops [29]. Because aerosol properties are unknown, 
it would be hard to accurately estimate the contribution to 
the adjacency effect caused by cloud-to-aerosol scattering, 
which can dominate at longer wavelengths. Thus extend- 
ing the shorter wavelengths correction into longer wave- 
lengths will improve estimates of the spectral signature of 
3D cloud enhancement. 

Next we revisit some features of the MODIS aerosol 
(MOD04) retrieval algorithm [20] that are relevant to this 
study. First, a cloud masking algorithm removes cloudy 
pixels within 20 pixel by 20 pixel boxes (with pixel size of 
500 m by 500 m each box is 10 km by 10 km). In addition 
to cloudy pixels, the darkest and brightest 25% are also 
discarded as likely affected by clouds [15]. The radiances of 
the remaining 500 m pixels are averaged and the average 
values are used later for multi-spectral aerosol retrievals. 
As a result, Eq. (2a) can be modified as 

R '3 d(2) = Rio(2) +4'(A), i = 1 , 2, . . . , N, (2b) 

where N is the number of remaining pixels in a 10 km by 
10 km box. Note that in order to be used for aerosol 
retrievals, a 20 by 20 pixel box must have at least 10 
(out of 400) remaining pixels, i.e. N > 10. Here we assumed 
that Rid does not vary for the selected N pixels. 

Let a and b be linear regression coefficients between 
R 3D values for two wavelengths in a single 10 km by 10 km 
box, i.e. 

R 1 30 ( 22 ) ^ GR^oOlij+b, i=l,2, (3a) 

We assume that Ai <A 2 (e.g., Ai=466 nm and A 2 = 855 nm). 
For the same two wavelengths, the linear regression for A's 
can be wr'tten as 

A\A 2 )^a A A\A^ + b A , 2 = 1,2, ...,N (3b) 

with regression coefficients a A and b A . Since it was 
assumed that Ri D is constant for the whole 10 km by 
10 km box, we expect 

a A ^ a. (4a) 

Indeed, 

4'(A 2 ) = R'3d(22)— Rid(2 2 ) ~ aR 1 3d(Ai) + b — Rid(2 2 ) 

= a[Rio(Ai ) +4'(Ai )] + b — Rid(2 2 ) 

= a4'(Ai) + aRiD(2i) — Rio(A 2 ) + b. 

Comparing this with (3b), we see that a A &a and 


b A ^ aRio(Ai)— RiD(A 2 )+b. 

(4b) 

Our hypothesis here is that 


b A ^ 0 

(5a) 

or, as follows from (5a) and (4b), 


^ 10 (^ 2 ) ~ Q^io(2i) + b. 

(5b) 
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In reality, because of aerosol small-scale spatial varia- 
bility (including changing of aerosol particles near clouds), 
aerosol (and surface) properties are not constant; thus Ru> 
varies inside a 10 km by 10 km box. (Note that for a special 
case of b= 0, Eqs. (3a) and (5b) are identical to Eq. (1).) 

The linear regression coefficients a and b (as well as a A 
and b A ) for given two wavelengths depend on many 
factors: solar and viewing geometry, aerosol spectral 
properties, cloud structure and even surface spectral 
reflectance. In general, a A ^ a and b A ^ 0. However, they 
can be equal approximately, and the validity (and accu- 
racy) of approximations (4a) and (5a) can be tested 
numerically with realistic models of broken cloud fields. 
Here we will use models of oceanic trade cumulus clouds 
and radiance enhancements in aerosol columns due to 3D 
cloud effects calculated with SHDOM [7]. 


As a result, if approximations (4a) and (5a) are accurate 
enough, the unknown 3D correction for longer wave- 
lengths can be approximated from the known correction 
for shorter wavelength as 

A app (A 2 ) = aA(M), (6a) 

and Ki D (A 2 ) can be approximated as 
R app \o(h) = R3D(h)~A a PP(X 2 ) = R 3 D(h)-aA(M), (6b) 

where a is determined from (3a). 

3. Data 

The UCLA LES model [22] was run for 26 simulations of 
marine boundary layer cumulus in a 20 km by 20 km 
domain with horizontal grid spacing of 62.5 m. The LES 



C 



Fig. 1 . (a) Distribution of aerosol height weighted with extinction, (b) Average cloud optical depth at 553 nm as a function of cloud fraction (with cloudy 
500 m pixels having optical depth 4 0.5) for the 26 LES scenes and (c) Number of clear pixels as a function of viewing zenith angles. 
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cloud microphysics assumes a fixed cloud droplet concen- 
tration (chosen randomly from an exponential distribution 
between 40 and 400 cm -3 ) and uses a two-moment 
drizzle/rain bulk microphysics scheme. Each simulation 
was initialized with a meteorological profile obtained from 
an ERA-lnterim reanalysis [6] column over the northeast 
Pacific between 2007-05-21 and 2007-09-23. The 3D 
liquid water content and relative humidity fields at 6 h 
simulation time were used. 

Forty representative aerosol profiles to use with the LES 
cloud fields were obtained from GEOS-5 [21] global fields 
over tropical and midlatitude oceans. GEOS-5 has five aerosol 
types: dust, sea salt, sulfate, black carbon, and organic carbon 
([4] and [5]), each with varying number of size bins. The 
complete aerosol optical properties were calculated from the 
GEOS-5 optical properties tables as a function of the LES 
relative humidity. The 40 aerosol profiles were chosen to 
have a wide range of 550 nm optical depth, single scattering 
albedo. Angstrom exponent, and extinction weighted mean 


altitude. The latter is illustrated in Fig. la. Each of the 40 
aerosol profiles was paired twice with one of the 26 LES 
cloud scenes to make a total of 80 cloud/aerosol scenes for 
this study. 

Radiances at the LES resolution were calculated with 
SHDOM [7] for these 80 scenes at the seven MODIS bands 
used for aerosol remote sensing (wavelengths of 0.466, 
0.553, 0.646, 0.855, 1.243, 1.632, and 2.119 urn). Rayleigh 
molecular scattering [2] and molecular absorption from 
water vapor and ozone were included up to the ~15.6 km 
domain top. An SHDOM ocean surface reflectance model 
with the LES surface wind speed was assumed. The solar- 
viewing geometry was obtained for the LES latitude, date, 
and the 13:30 local time of the Aqua overpass. SHDOM 
runs were made for 3D cloud/aerosol fields, 3D hydrated 
aerosol-only fields, and ID aerosol-only fields. The radi- 
ance differences between these runs show the effects of 
nearby clouds and variable humidity on the clear pixel 
aerosol radiances. MODIS 500 m pixels were simulated 




C 



Fig. 2. Probability density and cumulative functions of (a) a - a A , (b) b A and (c) b. 
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from the SHDOM radiances at 23 viewing directions (at 6° 
zenith angle spacing) appropriate for MODIS Aqua. Fig. lb 
shows the cloud fraction and average cloud optical depth 
for the 26 LES scenes. Finally, the MODIS Aerosol (MOD04) 
cloud masking procedure was implemented on the simu- 
lated MODIS pixels in 10 km blocks [15]. 

The total number of 500-m pixels that passed MOD04 
cloud masking for all 80 scenes is 100,188, while the total 
number of 10 km by 10 km boxes is 3154. Note that each 
scene consists of maximum four 10 km by 10 km boxes, 
but not all of them have more than ten 500 m by 500 m 
pixels that passed the cloud masking procedure, especially 
for high viewing zenith angles (Fig. lc). Hence the results 
provided below are slightly biased towards higher sun and 
lower viewing zenith angles. 

4. Results 

First we will test numerically the accuracy of the 
assumptions (4a) and (5a). Fig. 2a shows the probability 
density function (pdf) and the cumulative function of the 
differences between a and a A . It confirms that there is no 
bias (median=0.0) and in 80% of cases |3 - 1< 0.05. The 

average relative difference between a and a A is 2.4%. Fig. 2b 
illustrates the pdf of b A and its cumulative function. From 
the cumulative function we can see that in 94% cases p A \ 
<0.002, in 81% cases p A |< 0.001, and in 64% cases p A \ 
<0.0005. Finally, Fig. 2c confirms that the b-coefficient 
from Eq. (3a) is different from 0: in 78% cases p |4 0.01, in 
24% cases p |4 0.05 and in 5% p |4 0.1. 

Fig. 3a is a scatter plot of A app (X 2 ) vs. A(X 2 ) for 
A 2 =855nm. We assume here that the 3D correction for 
X\ is known with some uncertainties. We simulate them by 
multiplying A(X i) by (1 +k£), where k is the maximum 
uncertainty level and £ is a uniformly distributed random 
number between - 1 and 1. As follows from Eq. (6a), the 
multiplicative factor (1 +k£) applied to A(X i) can be also 
interpreted as the one applied to the regression coefficient 


a. But a is only an approximation to an unknown a A . Based 
on Fig. 2a, multiplication of a by (l+k£) does not bias it, 
increasing only slightly its standard deviation for 0<k 
< 0.4. As a result, the unbiased but randomly perturbed A 
( X \ ) does not affect much the accuracy of A app {X 2 y, Fig. 3a 
shows almost the same linear fit between A app (X 2 ) and A 
{X 2 ) in all three cases where (i) A(X i) is known exactly 
(k= 0), (ii) k= 0.2, and (iii) k= 0.4. 

In addition to scatter plot, the pdfs of A(X 2 )-A app (X 2 ) 
(1 +k%) and their cumulative counterparts are plotted in 
Fig. 3b. The cumulative functions for k=0 or k= 0.2 are 
very similar: quantitatively in these two cases 83% of all 
pixels have fi app (X 2 )-A(X 2 ){1 +k£) |< 0.001 and 63% pixels 
have ]fi app (X 2 ) — 4(2 2 )(1 + k£) |< 0.0005 (for k=0.4 these 
numbers are slightly lower: 81% and 59%, respectively.) 

The difference between R 3D (i 2 ) and its approximations 
(the corrected 3D and the ID ones) are plotted vs. R 3 d(X 2 ) in 
Fig. 4a. While the ID approximation is biased low, the 3D 
correction is uniformly distributed around the true values of 
R 3D (/l 2 )- This is confirmed with their pdfs plotted in Fig. 4b. 
Note that the pdfs of plotted in Fig. 4b are 

identical to pdfs of A{X 2 )—A app (X 2 ) plotted in Fig. 3b since, 

f^3D(^2) -K aPP 3D(>l2) = [Rm(X2)+A(X 2 )] - [Rw(X2) + A app (X 2 )] 
=A(X 2 )-A app (X 2 ). 

In addition to k=0, 0.2 and 0.4, Fig. 4b also shows the 
case with k= 0.8. As we see, the uncertainties in A(X 1 ) bias 
RWi(h) (or A app {X 2 )) only slightly. Quantitatively, means 
(and standard deviations) of the differences are 6.4 x 10 -5 
(0.0010), 6.4 x 10“ 5 (0.0011), 6.5 x 10" 5 (0.0012), 6.3 x 10“ 5 
(0.0015) for k= 0.0, 0.2, 0.4 and 0.8, respectively, while for the 
difference with ID approximation the mean is 0.0016 and the 
standard deviation is 0.0014. 

Finally, we note that as follows from Eq. (6b), our goal is 
to estimate the unknown # 10 (^ 2 ) (or A app (X 2 )) rather than 
K 3 DU 2 ). which is measured and does not need to be 
approximated. However, here the accuracy of R 1 FU 2 ) 
approximation was illustrated using the comparison 



Fig. 3. (a) Scatter plot of 4 app (2 2 ) determined from 4(2i) using Eq. (6a) vs. the "true" A(A 2 ) where Ai=466 nm and A 2 =855 nm. The unbiased uniformly 
distributed random uncertainties with maximum of 20% (blue) and 40% (black) have been added to A(Xf), and (b) Probability density and cumulative 
functions of 4(A 2 )— 4 app (A 2 ). 
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Fig. 4. (a) The difference between /?3 d(^ 2) and its two approximations: i?3lPU2)=RiDU2)+a4Ui) (red) and RioOte) (gray) vs. f?3D(/l2)- (b) Probability density 
functions of Rsd^)-^^-^) for 5 cases: A(Ai) known exactly; with maximum of 20%, 40%, 80% random unbiased uncertainties, and the case of RW’U 2 )= 
Rid(^ 2 )- As in the other figures, Ai=466 nm and ^2=855 nm. 


between Rfflifa) and R3DU2) which is the same as the 
difference between RiW(A 2 ) and K1DU2X 


5. Summary and discussion 

Ignoring 3D radiative cloud adjacency effects is one of 
the main sources of errors in remote sensing retrievals of 
aerosol properties in cloud-free columns [13,14,26]. Pre- 
venting these errors by excluding the retrieved aerosols 
near clouds will dramatically reduce the database and 
underestimate the aerosol radiative forcing, while includ- 
ing these areas may overestimate the forcing because of 
unaccounted cloud contamination. Many studies discussed 
this problem (e.g., [3,17,25]) but only few proposed ways 
to account for the effect [18,29,30,9-11]. So far none of 
them has reached a level of maturity for operational use. 

This paper combines the earlier approaches proposed in 
Marshak et al. [18] and Kassianov and Ovchinnikov [9]. First, 
based on Marshak et al.‘s paper, it is assumed here that we 
can successfully correct for 3D radiative effects at shorter 
wavelengths, where cloud-Rayleigh interactions dominate 
the radiation enhancement in cloud-free columns. Second, 
we use the Kassianov and Ovchinnikov approach, but instead 
of Eq. (1), we assume Eqs. (3a) and (5b). The early and new 
assumptions would be similar if b= 0, but it is not (see 
Fig. 2c). Third, in contrast to Kassianov and Ovchinnikov's 
approach, we do not retrieve AOT but only correct reflec- 
tances for the cloud adjacency effect. Thus we do not use the 
assumption of a power-law spectral representation of AOT 
that may be too restrictive. 

The proposed correction algorithm for a 10 km x 10 km 
box is the following: 

(i) At the shortest available wavelength Ai, get the correc- 
tion D(Ai) for cloud-Rayleigh interaction using the 
“two-layer” model of Marshak et al. [18]. 


(jj) For a longer wavelength A 2 , get a linear regression 
coefficient a relating R3DU2) to R3DU1) as in Eq. (3a) 
for all 500 m pixels kept by the MODIS operational 
cloud masking algorithm. 

(iii) Estimate A(A 2 ) as a product of a and A(A 1) and subtract 
it from measured R 3D {A 2 ) as in Eq. (6b). 

(iv) Repeat the above correction for all other wavelengths 
A4 A\ used in the MODIS spectral AOT retrieval 
algorithm. 

(v) Finally, apply the MODIS AOT retrieval algorithm to 
the estimated ID reflectances— that is, to the differ- 
ences between R 3D (/l) and A(A) for all bands. 


Note that the assumptions of dark surface and the 
aerosol layers below cloud top are applied only to our 
simple two-layer model (step (i)) and are not used with 
other steps of the algorithm proposed in this paper. The 
two-layer model has not been yet thoroughly validated. 
The results of validation will be reported elsewhere. 

To test our approach, we used 26 cumulus cloud fields 
(20 km by 20 km each) generated with the UCLA LES 
model, and combined them with 40 aerosol profiles based 
on GEOS-5 [21] global reanalysis. Radiances for 80 scenes 
were calculated using SHDOM [7] for two MODIS bands 
(Ai =0.466 Dm and A 2 = 0.855 Dm) and 23 viewing direc- 
tions. Radiances were averaged to 500 m resolution, and 
the MODIS aerosol cloud masking procedure was imple- 
mented on 10 km boxes [15]. The total number of 10 km by 
10 km boxes was 3154 and the total number of 500 m by 
500 m pixels that passed MOD04 cloud masking in all 80 
scenes and for 23 view directions was 100,188. 

The results showed that (i) for a variety of cumulus 
cloud scenes and aerosol profiles over dark ocean surfaces, 
the 3D correction due to cloud adjacency effect can be 
successfully extended from shorter to longer wavelengths 
and (ii) the 3D corrections at longer wavelengths are not 
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very sensitive to random unbiased uncertainties in the 3D 
corrections at shorter wavelengths. 
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